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ABSTRACT 

O ■ We investigate the evolution of a population of 100 dark matter satellites 

^ ! orbiting in the gravitational potential of a realistic model of M31. We find that 

^ ! after 10 Gyr, seven subhalos are completely disrupted by the tidal field of the host 

^ ! galaxy. The remaining satellites suffer heavy mass loss and overall, 75% of the 

mass initially in the subhalo system is tidally stripped. Not surprisingly, satellites 
with pericentric radius less than 30 kpc suffer the greatest stripping and leave a 
! complex structure of tails and streams of debris around the host galaxy. Assuming 

^ i that the most bound particles in each subhalo are kinematic tracers of stars, we 

find that the halo stellar population resulting from the tidal debris follows an 
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density profile at large radii. We construct B-band photometric maps of 
O ! stars coming from disrupted satellites and find conspicuous features similar both 

^ I in morphology and brightness to the observed Giant Stream around Andromeda. 

^ I An assumed star formation efficiency of 5-10% in the simulated satellite galaxies 

results in good agreement with the number of M31 satellites, the V-band surface 
O I brightness distribution, and the brightness of the Giant Stream. During the first 

I 5 Gyr, the bombardment of the satellites heats and thickens the disk by a small 

^ ' amount. At about 5 Gyr, satellite interactions induce the formation of a strong 

bar which, in turn, leads to a significant increase in the velocity dispersion of the 
^ ' disk. 

Subject headings: galaxies: interaction — methods: N-body simulations — cos- 
mology : miscellaneous 
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1. Introduction 

In the current cosmological paradigm, large-scale structure forms via hierarchical clus- 
tering wherein small systems composed of dark matter and gas merge to form larger objects. 
These systems originate from the primordial density fluctuation spectrum of cold dark mat- 
ter (e.g., Davis et al. 1985). The hierarchical clustering hypothesis leads to a picture in 
which subhalos are incorporated into larger systems. A number of processes such as tidal 
stripping and dynamical friction can lead to the destruction of substructure (White & Rees 
1978). However, many subsystems survive and remain in orbit about the parent halo. This 
phenomenon gives birth to a large population of small satellites orbiting around galaxies and 
clusters of galaxies (Moore et al. 1999; Klypin et al. 1999; Ghigna et al. 2000; Diemand et al. 
2004; Gao et al. 2004; Benson 2005). Typically this population consists of several hundred 
subhalos and contributes about 10% of the total dark halo mass in a galaxy (Font et al. 
2001; Ghigna et al. 2000). 

From an observational point of view, there is no clear evidence that such a large popu- 
lation of dark matter satellites exists. There are only ^ 40 observed satellites in the Local 
Group out of which 13 belong to the Milky Way (Mateo 1998). This number is an order 
of magnitude less than predicted from cosmological simulations because the stellar content 
and hence surface brightness is below current detection limits. How then is it possible to 
reconcile the low number of observed satellites with models predictions? There are several 
possibilities that we examine below. 

The satellites may still be there but difficult to observe. Many of the satellites may 
be associated with high-velocity clouds (Blitz et al. 1999) or may not be detected because 
the stellar component is not important enough and the surface brightness of these objects is 
below current detection limits. There are also theoretical reasons why star formation may 
be suppressed or inefficient in low mass satellites. Several studies explore mechanisms that 
should operate in the early stages of galaxy formation to suppress star formation in low- 
mass satellites. Gas ejection by supernovae-driven winds from the shallow potential wells of 
satellites may quench star formation (e.g. Dekel & Silk 1986). Also a strong intergalactic 
ionizing background can prevent the collapse of gas into low-mass systems (Thoul & Weinberg 
1996) with circular velocities Vc ^ 30 km as revealed in numerical simulations. Recently, 
Minchin et al. (2005) claim to have detected an HI source in the Virgo Cluster that could 
be a dark halo without an accompanying bright stellar galaxy. If the observed system is 
bound, this implies a dark halo mass up to 10^^ Mq. In addition, the census of the Milky 
Way satellites may be incomplete at low latitude due to obscuration. Willman et al. (2004) 
estimated a 33% incompleteness in the total number of dwarfs due to obscuration which 
represents only a small increase of the total number compared to cosmological predictions. 
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AU these possible explanations and considerations are plausible solutions to the sateUite 
abundance problem. 

A large sateUite population may produce a strong dynamical effect and modify the 
structure and kinematics of the galactic disk. Many studies look at the dynamical response 
of the host galaxy from an infaUing sateUite (Toth & Ostriker 1992; Quinn et al. 1993; Walker 

et al. 1996; Huang & Carlberg 1997; Velazquez & White 1999; Benson et al. 2004). Mihos 
et al. (1995) show that the accretion of a low- mass satellite by a disk galaxy can generate a 
strong bar. This bar can buckle vertically and eject disk material out of the disk giving rise 
to an X-structure similar to structures observed in some bulges (e.g. Whitmorc & Bell 1988). 
Angular momentum transfer from the satellite to the disk may also occur during the infall 
phase. Huang & Carlberg (1997) demonstrate that sateUites with a mass 0.3Mdisk change 
the orientation of the disk up to 10° and generate warps as the satellites, under dynamical 
friction, sink toward the center of the galaxy. In addition, Quinn & Goodman (1986) and 
Quinn et al. (1993) note that the effects on the vertical structure of the disk are not uniform 
across the disk. They observe some flaring of the disk at large radii, i.e. particles at large 
radius are orbiting at larger z than those at small radius. As suggested by these authors, a 
thick disk produced by minor mergers should have a scale height that increases with radius. 
Repeated interactions of a large population of satellites with a stellar disk also lead to disk 
heating. Quinn & Goodman (1986) shows that the disk is not heated isotropicaUy by the 
infall of a sateUite. They note that the heating is largest near the center of the disk and 
that most of the kinetic energy is distributed in the disk plane, i.e. receives a small 
fraction of the available energy. In fact, most of the radial and azimuthal heating comes 
from the disk spiral response to an infaUing satellite (Quinn et al. 1993). In addition, the 
disk spreads radially due to an angular momentum exchange between the satellite and the 
disk. Velazquez & White (1999) demonstrates that the heating and thickening rate of the 
disk differs for satellites on progradc and retrograde orbits. The former heat the disk and the 
second tilt it. Benson et al. (2004) show that there are significant differences between heating 
rates for prograde and retrograde orbits and those differences are amplified by increasing the 
mass and concentration of the sateUites. As noted by some authors (Benson et al. 2004; 
Font et al. 2001), only satellites with orbits that bring them close to the center of the galaxy 
affect the disk significantly. 

Most of the results described above were based on controlled numerical experiments 
where a single satellite interacts with a "five" disk. The exception was the pioneering study 
by Font et al. (2001) who examine the evolution of a stellar disk embedded in a system of 
several hundred subhalos. The work presented here improves upon that study in a number of 
ways. First, the mass resolution of our simulations is a factor of 250 higher than in the Font 
et al. (2001) study. This is high enough to adequately follow the evolution of the satellites 
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and disk and to suppress two-body effects. Second, tfie base model for the galaxy is the self- 
consistent equilibrium model of M31 designed to fit the observed photometric and kinematic 
data for the Galaxy. Finally, the assumed mass distribution of the subhalo population 
and internal structure of individual subhalos are motivated by results from cosmological 
simulations. 

Om aim is to study the effect of tidal stripping on the satellite population. In particular, 
we attempt to quantify the number of satellites that survive with a detectable stellar surface 
brightness. (Our subhalos have a single component, but we can infer the surface brightness 
by modelling the stellar content, as described below.) In addition, we study the stellar halo 
formed from tidal debris paying particular attention to extended structures. As we will see, 
these structures are similar to observed features of M31 such as the Giant Stream (Ibata et 
al. 2001). 

In §2, we describe our N-body model of M31 and describes a control experiment with- 
out sateUites designed to test the model's inherent stability. In §3, we describe the initial 
conditions of our satellite system. In §4, we discuss the evolution of the sateUite population 
and follow with a discussion of the hypothesis that the metal-poor stellar halos arise in the 
tidal stripping of satellite galaxies. Our main conclusions are summarized in §5. 

2. An N-body model of M31 

For our experiments, we use a self-consistent numerical model of M31 derived from a 
composite distribution function (Widrow & Dubinski 2005) (hereafter the WD models). In 
general, the WD models are axisymmetric equilibrium solutions of the coUisionless Boltz- 
mann equation which may be subject to non-axisymmmetric instabilities. We choose their 
model M31a which provides a good match to the observational data for M31. Numerical 
simulations of this model assuming a smooth halo find that the system is stable against the 
formation of a bar for 10 Gyr. Therefore any heating or instabilities that are observed in 
simulations where part of the smooth halo is replaced by a subhalo population should be the 
result of the interactions between the disk and the satellite system. 

The model consists of an exponential disk, a Hernquist-model bulge and a truncated 
NFW halo with plausible mass-to-light ratios for the disk and bulge and simultaneous fits 
to the surface brightness profile, rotation curve and bulge velocity dispersion and rotation 
with an NFW halo that extends to the expected virial radius. These models contrast with 
other studies of infaUing sateUites (e.g Benson et al. 2004; Velazquez & White 1999) that use 
disk galaxy models that must be pre-relaxed from an approximate equilibrium initial state 
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(Hernquist 1993). The M31 model has the main advantage that the initial state is formally 
an equilibrium solution to the collisionless Boltzmann equation since it is derived from a 
distribution function constructed from integrals of motion. 

Table 1 contains the physical parameters of our Andromeda model. The bulge follows 
the Hernquist profile 

pH(r) = 3 (1) 

{r/n) (1 + r/n) 

while the halo is modeled according to the NFW density profile (Navarro et al. 1996) 

Pnfw ~ ~ 2 ■ (^) 

r/rs (1 + r/rj 

The disk distribution function is taken from Kuijken & Dubinski (1995). In brief, the surface 
density profile follows an exponential with scale radius Rd and truncation radius Rout- The 
vertical component is given by secli' {z / 2zd) where Zd is defined as the scale height. 



2.1. Control Experiment Results 

We first run a control experiment of the M31 system for 5 Gyr using a parallel trcccode 
(Dubinski 1996). The softening length is 15 pc. We use 10,000 equal time steps and the 
total energy is conserved to within 0.7%. We take the 4.25 Gyr snapshot of the 35M- 
particle control model as the galaxy initial conditions for our run with satellites. Most 
of the initial transient spiral instabilities have died away by that time and the increase in 
velocity dispersion after that point is very small. Starting at this initial time minimizes 
the contamination of disk heating by initial transient spirals. We also run other control 
experiments at lower resolution (350K and 3.5M particles) to look for convergence. The 
simulations are run using CITA's McKenzie cluster (Dubinski et al. 2003). 

After taking the 4.25 Gyr snapshots as initial conditions, we let the galaxy evolve for 
10 Gyr. Figure 1 shows the evolution of the disk velocity dispersion for the 35M model. 
Outside 5 kpc, (7^ and increase by 5-10 km s~^ while cr^ is virtually unchanged. Inside 
5 kpc, all components increase significantly. In fact, most of the heating occurs in the first 
few billion years where spiral patterns arising from swing-amplified noise heat the radial and 
azimuthal components of the velocity ellipsoid. Nevertheless, the heating is quite small and 
negligible in the z direction. The disk scaleheight evolution for the same model is shown in 
Figure 2 and found to be minimal in the inner region of the disk. There is some flaring on 
the edges, but the increase is only around 100 pc. The velocity dispersion ellipsoid ((Tr,(T</,,(72) 
at an equivalent solar radius ol R — 2.5Rd changes from (17.8, 24.6, 18.9 ) km s~^ to (20.5, 
29.0, 19.1) km after 2.5 Gyr. In comparison. Font et al. (2001) let their disk relax for 
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3.5 Gyr before introducing the satellite population and their velocity dispersion ellipsoid 
changes from (31, 27, 26) km to (43, 30, 33) km s~^. Because we use a large number of 
particles for the galaxy (20M - halo, lOM - disk, 5M bulge), the mass of the halo particles 
are quite small and artificial disk heating due to halo particles bombarding the disk is almost 
negligible. 



The properties of the initial satellite population arc motivated from the latest results of 
ACDM simulations of cluster and galaxy formation. Since the properties and distribution 
of subhalos do not differ significantly on galaxy and cluster size scale (Moore ct al. (1999), 
L. Gao private communication), the different relations at cluster size scale (number density 
profile, slope of the mass function, total mass ratio, etc.) can be rescaled down to the galactic 
scale. 



We use the following mass function in agreement with recent AGDM simulations by 
(Gao et al. 2004) : 



Here M represents the number of sateUites per unit parent mass. We distribute the satellites 
between Mgat/Mhaio = 1-5 x 10~^ to Mgat/Mhaio = 0.02 in 15 different logarithmic mass bins. 
With Mhaio = 5.8 X 10^^ M©, we get n ~ 28 and a mass fraction of 0.024. In Gao et al. 
(2004), Figure 12 suggests that most of the satellites for a M31-sized galaxy are in place 
between 2; = 0.7 and 1.1. Their Figure 13 shows that satellites accreted aX z = 1 lose of 
order 70-80% of their mass by z = 0. With these results in mind, we start our simulation 
with a total mass fraction of 0.1 and 100 satellites. The sample of 100 satellites is large 
enough to provide the correct stochastic treatment of disk heating and will provide a smooth 
distribution of tidal debris. Figure 3 shows the comparison between the expected satellite 
cumulative mass function, the integral of Equation (3) normalized to 0.1 for 1.5 x 10~^ 
< -^satZ-^haio < 0.02, and the data points generated for our population composed of 100 
satellites. 



3. Initial conditions of the satellite population 



3.1. 



The mass distribution function 




(3) 
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3.2. Spatial Distribution of Satellites 

We assume that the spatial distribution of satelUtes in the parent halo is spherically 
symmetric. Recent cosmological simulations (e.g Diemand et al. 2004; Gao et al. 2004) 
suggest that the number density profile of satellites is less centrally concentrated than the 
dark matter halo. Gao et al. (2004) found that the cumulative number density profile of 
subhalos in a Milky Way size halo is well fitted by 

Ar/ N Ar 2 75 (1 + 0.244a;,) 

where x = r/r20Q, Xs = rs/r2oo and Nfot is the total number of satellites within r2oo- We 
sample that cumulative number density profile according to this formula with the further 
assumption that the spatial distribution of satellites is independent of the satellite mass. 



3.3. Satellite orbital distribution 

According to cosmological simulations, satellite orbits are isotropic in the center and 
slightly radially biased in outer regions (e.g. Benson (2005); Diemand et al. (2004)). The 
anisotropy parameter 

varies with radius as P{r/r 200) ~ 0.35r/r2oo for the satellite population (Diemand et al. 2004). 
We use a constant value oi P — 0.3 over the range < r/r2oo < 1 as a good approximation. 
The velocity components ((7^, ae, a^) are determined using the Jeans equation (Binney & 
Tremaine 1987) 

where n is the satellite number density profile and $ is the gravitational potential. Here 
we assume that the gravitational potential of the host galaxy can be approximated by an 
NFW density profile and the gravitational potential generated by the sateUite population 
is a small perturbation of the host gravitational potential, i.e. we neglect the contribution 
of the satellites. The NFW approximation is vahd everywhere, except within ~ 30 kpc 
where the bulge and disk modify the potential significantly (r < 30 kpc) . We solve equation 
(6) numerically for o"r(r) and assign the velocity {vr,V0,V(j)) at position r assuming they are 
random variates with dispersions {ar, Pctq, Pcr^). 

To validate the equifibrium of the sateUite distribution within the M31 model galaxy, 
we replace the satellites by test particles and compute the forces assuming a rigid potential. 
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So long as the satellites are test particles neither dynamical friction nor tidal stripping come 
into play. In Figure 4, we show the evolution of the cumulative number density function 
of satellites. There is an evolution of the distribution, but this may be attributed to small 
number of satellites and the inexactness of using the Jeans equations dispersion profiles to 
approximate the distribution function. Despite this slight evolution, we believe that our 
initial distribution could represent a population in quasi-equilibrium with the host galaxy 
and is therefore suitable for our experiments. 



3.4. The internal structure of the satellites 

Individual satellites are modeled as single cuspy NFW halos truncated by the tidal 
field of the galaxy model. The WD model for an isolated NFW halo has an adjustable 
parameter 7 that can introduce an arbitrary tidal radius akin to the the tidal radius of a 
King model. Satellite models can therefore be created that have an r"^ cusp with NFW 
behavior at the extremities and a self-consistent cut-off. Many previous studies used the 
King model for the sateUites (e.g. Benson et al. 2004; Velazquez & White 1999; Huang 
& Carlberg 1997) in sateUite-disk interactions. Since NFW sateUites are cuspy, we might 
expect them to be more robust to tidal interactions than King-model satellites and therefore 
produce more damage to the disk. Wc use 15 different satellite models to cover the mass 
range 1.5 X IQ-^ < Msat/A4aio < 0.02 logarithmically This mass range corresponds to 
8.75 X 10^ < Ms^i/Mq < 1.17 X 10^°. As a comparison, the Sagittarius dwarf and Large 
Magellanic Cloud, our galactic companions, have estimated masses of (2-5) x 10^ Mq (Law 
et al. 2005; Ibata et al. 1995) and 2 x 10^° Mq (Schommer et al. 1992). Note that the mass 
of the LMC corresponds to the upper hmit of the mass function. M31's satellites M32 and 
NGC205 have masses of 2.1 x 10^ and 7.4 x 10^ Mq respectively (Mateo 1998). Our mass 
range includes most of the observed satellites in the Local Group. 

For simplicity, the truncation of each satellite is determined by the length of the tidal 
radius at the mean apocentric radius of the satellite system (r = 50 kpc). Here we assume 
that the tidal radius of the sateUites r^, is defined by the Jacobi approximation (Binney & 
Tremaine 1987). The radial extent of the satellites is determined by comparing the mean 
density of the satellite inside the tidal radius and the mean density of the halo inside r = 50 
kpc : 

P~s{rt) ~ 3p^(r = 50kpc) (7) 

where 'Jhi'f^t) is the satellite mean density inside the tidal radius and 'Phi'r = 50kpc) is the 
halo mean density inside 50 kpc. Satellites within r = 50 kpc will over fill their tidal radii 
while those beyond will lie well within it. We determine that satellites initially located inside 
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50 kpc will lose a cumulative 5.3 x 10^ Mq due to the tidal radius overestimate. This accounts 
for an "artificial" tidal stripping at the beginning of the simulation and corresponds to less 
than 1% of the total satellite mass. A true satellite system is never really in equilibrium since 
dynamical friction leads to continual tidal erosion so our choice of a mean satellite orbital 
radius for estimating the tidal radius is a compromise. We see below that there is indeed 
a transient startup where satellites that overfill their tidal radii are quickly stripped but by 
an amount that is small (less that 1%) compared the total stripping over the course of the 
satellite system evolution. 

Table 2 contains the main characteristics of each satellite model. Wc model our satellites 
assuming a shallow power law distribution for the concentration of the satellites (Navarro et 
al. 1996, 1997). The internal properties of the satellite halos found in cosmological simula- 
tions are not well determined, especially for the low-mass end where the concentration and 
density profile are not well constrained because of the poor mass resolution. 

3.5. Evolution of individual satellites 

To compute the mass, center-of-mass position and density profile evolution of each 
infalling sateUite, we use the technique described in Benson et al. (2004). This algorithm 
identifies the particles that are bound to the sateUites, computes the center-of-mass of the 
system, and iterates until the total mass converges. Typically, the criterion we use for 
convergence varies by less than 0.5% in mass and most of the steps require only a few 
iterations. The method is straightforward and offers a good alternative to the friends-of- 
friends algorithm which is not appropriate when the satellites are numerous and characterized 
by different scale lengths. 

4. Results 

We perform two simulations at low and high resolution to test for numerical convergence. 
In the low resolution run, each of the one hundred satellites is represented by 1000 particles 
and put in orbit around an M31 galaxy model with lOOK disk particles, 50K bulge particles 
and 200K halo particles. In the high resolution simulation, particle numbers are increased 
by a factor of 100. The population of satellites for both simulations have the same mass 
function, spatial distribution and orbital velocities. Both simulations are run for 10 Gyr 
(20,000 equal time steps) using a parallehzed treecode (Dubinski 1996) with opening angle 
parameter 9 — 0.9 (quadrupole order) and e = 25 pc (e = 15 pc) for the 450K (45M) run. 
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For the 45M-particle run, the energy is conserved to within 0.4 %. Unless otherwise stated, 
all the results presented in this section will be for the high resolution version. 

4.1. GalcLxy Evolution in Three Acts 

4.1.1. Act I. The First Orbit 

During the first few billion years, the satellites accomplish their first complete orbits 
and leave behind an interwoven web of tidal streams. Some of these streams extend beyond 
the virial radius of the host galaxy. One of the most interesting features of the first two 
billion years is the presence of shell structures with clearly defined edges similar to those 
seen in some elliptical galaxies. These shells are produced via phase wrapping (Quinn 1984; 
Hernquist & Quinn 1987) and tend to disappear quickly due to phase mixing. During this 
period, the disk remains quiet with only a small amount of heating. 

4.1.2. Act 11. Bar Formation 

At about 5 Gyr, a stunning phenomenon appears in the disk : the unexpected formation 
of a bar. We can associate its creation with the interaction of the disk and the satellite system 
since no bar formed in the control run after 10 Gyr. The strong bar is responsible for most 
of the disk heating after 5 Gyr and creates a puffy vertical structure. 

4.1.3. Act III. Anticlimax : The Quiet End 

Following bar formation, the disk and satellite system evolution is more gradual. After 
two or three pericenter passages, the satellites are stripped significantly. The galactic halo 
is enriched with tidal debris and the initial streams become more mixed. The final state 
of the debris is that of a spheroid of about the size of the Galactic stellar halo. During 
that time, the bar suffers a bending instabihty that gives rise to the buckhng mode (Raha 
et al. 1991). This process generates a conspicuous X-design in the disk when viewed edge- 
on and may point to the origin of the peanut-shaped bulges observed in many galaxies 
(Bureau & Freeman 1999; Kuijkcn & Mcrrificld 1995; Whitmore & Bell 1988). During the 
final two billion years, the bar remains the dominant feature of the disk but slowly spins 
down. By the end of the simulation, a typical satellite completes five to seven orbits and 
loses a significant fraction of its initial mass. Snapshots of the sateUite population and the 
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host galaxy are shown in Figure 6 (see also the animation of the disk galaxy evolution at 
http://www.cita.utoronto.ca/~jgauthier/m31). 

In Figures 7 and 8 we show the evolution of the disk velocity dispersion ellipsoid and 
scaleheight for the high-resolution run. As noted in Figure 7, the formation of the bar 
around 5 Gyr produces a sudden jump in velocity dispersion for disk particles. The same 
event triggers the scaleheight increase in Figure 8. Further details of the disk evolution will 
be discussed in a forthcoming paper. 



4.2. Evolution of the Satellite System 

The satellites are initially distributed isotropically in space between approximately 50 
kpc to 250 kpc (Ri r2oo)- Although most of them survive the strong tidal interactions in the 
vicinity of the disk, they create extended tidal streams and shell structures. These streams 
and shells are particularly obvious in the first few billion years, but tend to lose their sharp 
edges as phase-mixing proceeds. In Figure 9, we show the evolution of the two most massive 
and three least massive satellites. The plots, which show the evolution of the satellites 
position and density profile demonstrate that dynamical friction has little effect in bringing 
the satellites close to the center of the galaxy. Figure 10 shows that there is no significant 
change to the number density profile of satellites after 10 Gyr. Because of the very steep 
slope of the mass function, most of our satellites have a mass of 1.5 x 10~^ Mhaio and do 
not feel the effects of dynamical friction. In fact, plots of the satellite orbital decay indicate 
that the value of the apocenter radius decreases very slowly for almost all satellites, except 
for the very massive ones (see Figure 9). 

The inner slope of the satellite density profile remains the same over the simulation (see 
Figure 11). As noted by Hayashi et al. (2003), it is possible to describe the structure of a 
stripped halo by modifying the NFW profile : 

Pir) = rK^Arys^NFw (8) 

where ft is interpreted as a reduction in the central density of the profile and rte is an 
effective tidal radius that describes the cutoff due to tidal forces. Comparisons are hard to 
make with the work done by Hayashi et al. (2003) because our sateUites suffer an initial 
truncation to avoid a divergent mass profile. In some way, our initial satellites density 
profiles look quite similar to the final profiles of their subhalos. Nevertheless, we show, in 
Figure 11 the best fit of equation (8) for typical profiles of three different satellites mass bins. 
Equation (8) provides a good fit of the final profile especially for the low-mass satellites with 
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Msat/Mhaio < 5.5 X 10 ^. For the more massive ones, the fit tends to overestimate the mass 
loss at large radii. 

Though only seven satellites are completely destroyed by the end of the simulation 
most of the remaining satellites are stripped significantly. As shown in Figure 12, the time 
dependence of the total mass bound in satellites can be described by two distinct phases. 
During the first 4 Gyr, the satellites lose approximately half their mass. The time dependence 
of the mass in the system is well-fit by an exponential decay : 

M,at oc e-°-^^='*/*V2 (9) 

where ti/2 = 3.5 Gyr. The second phase, from 4 to 10 Gyr, is quiet with ti/2 = 9-3 Gyr. 
During their first complete orbit, satellites are severely stripped and lose their outer mass 
layers as their size becomes limited by the tidal radius at the pericentric passage. For the 
last 6 Gyr, the satellites radial extent is well constrained by the tidal field of Andromeda and 
a smaller mass loss occurs at each pericenter passage. Glearly these numbers are affected by 
our initial conditions and especially by the constraints on the satellites size. The position 
at which we compute the tidal radius of our satellites (50 kpc) is somewhat arbitrary and a 
different value could lead to different results. As shown in Figure 5, the distribution of the 
pericenter radii peaks at 50 kpc or so. Gomputing the tidal radius of the satellites at this 
position gives a good approximation of the radial extent of satellites who have completed 
several orbits around the galaxy (i.e. suppose to be in equilibrium with the host galaxy). 
Nevertheless, decreasing the radius to 30 - 40 kpc means increasing the concentration of the 
satellites to unrealistic values. It appears that the tidal field of the underlying halo potential 
is very strong and satellites faUing inside 30-40 kpc must have a very high concentration to 
survive multiple pericenter passages. On the other hand, computing the tidal radius at say, 
100 kpc and constraining the satellites to have a size less than rtidai means that the satellites 
will reach the disk galaxy with most of their mass stripped. We think that a value of 50 kpc 
(4 Vs) is a reasonable choice. 



4.3. Absence of Holmberg Effect 

Holmberg (1969) showed that there is a tendency for satellite galaxies to congregate near 
the poles of the spiral host galaxy. His observations were then confirmed by Zaritsky et al. 
(1997) who showed that sateUites located at large projected radii of isolated disk galaxies are 
aligned preferentially along the disk minor axis. However, recent observations by Brainerd 
(2005) on a sample of SDSS galaxies showed that satellites are preferentially aligned with 
the major axis of the galaxy. This result contradicts Holmberg's previous observations. In 
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this paper, we examine if either of these conclusions are detected in our sample of evolved 
satellites and if a strong dynamical interaction between the disk and satellites at low latitudes 
could explain the Holmberg Effect. If dynamical friction is more important at lower latitudes, 
one would expect that satellites on almost coplanar orbits to sink on shorter timescales than 
satellites on polar orbits. That would lead to a deficit of satellites at low latitudes. In 
Figure 13, we show the distribution of satellites as a function of cos 6*. The results are 
consistent with a uniform distribution in cos^, that is a spatially isotropic distribution. In 
other words, we do not detect any Holmberg or anti-Holmberg effect in our simulation. We 
conclude that the anisotropic distribution of satellites observed around galaxies does not 
come from a dynamical interaction with the host galaxy but probably originates from the 
galaxy formation initial conditions. 



4.4. Outer halo stellar density profile 

An important aspect of the satellite galaxy problem is their detectability. Our interest 
is in finding the position of the stars coming from disrupted satellites and making predictions 
about their distribution and luminosity. We assume that a good proxy for the initial position 
of stars are particles deep in the potential well of the host satellite. For our purpose, we 
assume that the 10% most bound particles are kinematic tracers of the stellar population in 
each satellite. The rationale for it is simple : as gas cools down, it loses energy and sinks 
in the potential well of the host satellite. Thus, we expect that most of the stars have large 
binding energy. Once these particles are identified, they are labeled and followed during the 
simulation. Each point particle represents a population of stars. To assign these particles 
a stellar mass, we assume a baryonic mass fraction and a star formation efficiency. We 
normalize the mass of these particles so that they correspond to a baryonic mass fraction of 
0.171 (Steidel et al. 2003) and a star formation efficiency, the fraction of baryonic mass turned 
into stars, between 1 and 10% (Ricotti & Gnedin 2005). We assume a simple scenario in 
which the baryonic mass fraction and star formation efficiency are the same for all satellites. 

We show in Figure 14 the spherically averaged density profile of "stars" . We include the 
stars that have been tidally disrupted from their host satellites and the ones that have not 
been displaced from the center of the potential well of individual satellites. The contribution 
from satellite stars dominate for r/r2oo > 0.2 and the density profile of stellar material at 
larger radii can be well fit by an r~^'^ power law that agrees with globular clusters system 
(Harris 1976) of our Galaxy. This profile is also in good agreement with observations of the 
Milky Way's metal-poor stellar halo (Morrison et al. 2000; Chiba & Beers 2000; Yanny 2000) 
and with recent N-body simulations by Abadi et al. (2005) and Bullock & Johnston (2005) 
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who show that the density profile of the accreted stars goes as p oc r~°, a ~ 3-4. In addition, 
recent analysis of 1047 SDSS edge-on disk galaxies by Zibetti et al. (2004) demonstrate the 
presence of stellar halos with spatial distribution that is well described by a power law 
p oc r~^. Our simulation predicts that, after 10 Gyr, the satellites contribute a total stellar 
mass of 1.2 x 10^ Mq. About 20% of this mass is found inside the edge radius of the disk 
and represents less than 1% of the total stellar mass of the disk and bulge components. The 
overall contribution of the satellite debris to the total galactic stellar mass is about 1%. 



4.5. B-band photometric maps : Streams 

In order to convert surface density (which one gets from the N-body distribution) to 
surface brightness, one requires the mass-to-light ratio of the stellar population. This in 
turn requires a stellar evolution model and we use those by Bruzual & Chariot (2003)^. 
The resulting mass-to-light ratios in the B-band at different epochs along with the baryon 
fraction and the star formation history are listed in Table 3. We assume a constant T = 7.6 
B-band mass-to-light ratio for M31 (Berman 2001; Faber & Gallagher 1979). 

We generate 2000 x 2000 pixels photometric maps showing a 10° x 10° field centered 
on M31 so that each pixel has a corresponding plate-scale of 18 arcsec. To reduce shot noise, 
we smooth out our "point mass stellar populations" with a 5 pixels (1.5 arcmin) gaussian 
window. The final results are presented in Figure 17 and represent B-band photometric 
maps of the stellar population expected from disrupted satellites galaxies taken at different 
times. We rotate our N-body model to fit the orientation of M31 on the sky. For each panel, 
we vary the star formation efficiency and produce three different maps. 

We now examine the tidal streams in our simulations. Figure 18 presents a close- 
up of photometric maps at 3.5 and 5.5 Gyr. At 3.5 Gyr, there is an obvious "bridge" of 
stars connecting two subhalos and M31 in the upper map. This feature is similar both in 
morphology and brightness to the M31 giant stream detected in a number of surface density 
maps (Ferguson ct al. 2005, 2002; Ibata et al. 2001) and constitutes the brightest stream 
detected in our photometric maps. The mean surface brightness of our giant stream proxy 
is about 28.5 mag arcsec"^ (in the B-band) and has comparable brightness to the actual 
measured value of //^ 30 ± 0.5 mag arcsec"^ (Ibata et al. 2001). Even though we generate 
maps in the B-band, our photometric data in the B-band are very similar to those one would 
obtained in the V-band because, for a few billion year old simple stellar population, the mass- 
to- light ratio in the B-band is similar to the one in the V-band (//& ~ pLy in the simulations) . 



See http : / /www. cida. ve / "bruzual /bc2003 
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Other streams have also been observed around M31 and have similar brightnesses. Zucker 
et al. (2004) find a 3° overdensity of luminous red giant stars (Andromeda NE) having a 
central g-band surface brightness of 29 mag arcsec"^. In addition, McConnachie et al. (2004) 
observe an arc-like overdensity of blue, red giant stars in the west quadrant of M31. This 
tail has a fiy = 28.5 ± 0.5 mag arcsec^^. Similarly, the Gl Clump, a stellar overdensity that 
is located 30 kpc along the southwestern major axis could be associated by an overdense 
clump locate at the lower right edge of Andromeda in bottom picture of Figure 18. 

4.6. Surface Brightness Profile 

In Figure 16 we show the surface brightness profiles of the stars no longer bound to the 
satellites and contributing to the metal-poor halo of M31. The plots arc drawn assuming a 
star formation efficiency of 10% and a simple stellar population formed 1 Gyr prior to the 
beginning of the simulation. Comparisons with actual observations of the surface brightness 
profile of M31 are quite challenging because the extremely low surface brightnesses involved 
pose a significant problem for observers (typically, 7-8 mag higher than the sky). One of the 
deepest surveys of M31 halo surface brightness was carried out by Irwin et al. (2005). These 
authors have been able to go down to yUy ^ 32 mag arcsec"^ at a projected radius (along 
the minor axis) of 4° (55 kpc). At that radius, they obtain a V-band surface brightness of 
about 30-31 mag arcsec^^ which is close to our B-band estimate of 29-30 mag arcsec^^. In 
Figure 16, a star formation efficiency of 5% would decrease the surface brightness at all radii 
by about -1-0.75 mag arcsec"^ and leads to a better agreement with the value measured at 
55 kpc. 

The value of the surface brightness for R < 50 kpc flattens at about 28 mag arcsec^^. 
Dynamical friction only weakly affects the trajectories of the satellites and thus, little mass 
is being deposited in the center of the galaxy. Other studies (e.g. Abadi et al. (2005)) show 
a constant increase of the surface brightness profile down to a radius of about 10 kpc. In 
many simulations, the satellites sink in the center of the galaxy and get heavily stripped 
within 30 kpc and deposit large amounts of mass in the center. It is not the case in this 
simulation. 



4.7. Satellites detected around M31 

A way of constraining the star formation history of the satellites is by counting the 
number of sateUites detected in the close vicinity of M31 as a function of the star formation 
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efficiency and epoch of birth and compare their number and properties with the current 16 

M31 satellites claimed by McConnachie & Irwin (2006). The faintest of these satellites is 
Andromeda IX {nv = 26.8 mag arcsec"^) discovered by Zucker et al. (2004). The survey is 
complete for satellites with fiy < 26.8 mag arcsec^^. According to the results listed in Table 
4, it appears that a star formation rate of about 10% combined with an epoch of formation 
between 3.5 and 1 Gyr prior to the beginning of the simulation can account for the actual 
number of sateUites detected around M31. For the 5.5 Gyr snapshots, a star formation rate 
of 10% results in exactly 15 detected sateUites with //y < 26.8 mag arcsec"^. On the other 
hand, a SFR of 1% is insufficient since it reduces to the number of sateUites detected below 
the 26.9 mag arcsec"^ hmit. 

We expect the most massive satellites to have a light distribution similar to M32 or NGC 
205. The galaxies have, respectively, a central B-band surface brightness of approximately 
17.5 mag arcsec"^ and 20 mag arcsec"^ (Choi et al. 2002) . The brightest simulated satellites 
located in the vicinity of M31 have a central B-band surface brightness of about 23 mag 
arcsec"^ assuming a simple stellar population formed 1.5 Gyr before the beginning of the 
simulation and and a star formation efficiency of 10%. The fact that our simulations are 
dissipationless is the cause of this difference. In a more realistic case, gas cools down and sink 
deep in the potential well increasing the central density of stars. Our simulations provide 
lower limits to the value of the central surface brightness one would expect. However, a 
central value of 22-23 mag arcsec"^ is consistent with a recent study of dwarf galaxies in the 
Fornax cluster (PhiUipps et al. 2001) and the Andromeda (I- VII) dwarfs series (McConnachie 
& Irwin 2006). These authors find that the majority of the dwarfs in their sample have a 
surface brightness spanning 20 - 24 mag arcsec~^ (see their Figure 1). Many of our satellites 
would fall in that range of central brightness. 

The surface brightness at the half-light is a more representative measurement of the over- 
all dwarf light distribution and does not suffer much from dissipation effects. Our brightest 
satellites have A*e//,b — 25.3 mag arcsec"^. By comparison, M32 has a half-light surface 
brightness fJ'eff,B ~ 20 mag arcsec"^ (Choi et al. 2002). A more recent episode of star 
formation or a continuous star forming period could explain that difference since younger 
stars would increase the B-band ffux emission. Our simple model assume that all stars have 
the same age and in this case, formed 1.5 Gyr before the beginning of the simulation. We 
conclude that our brightest satellites show photometric differences with the most massive 
companions of M31. However, our simulations are purely coUisionless and including dissi- 
pative processes would increase the central and half-light radius luminosities in such a way 
that the brightest simulated sateUites would be similar to M32 and NGC205. 

Our choice of a satellite mass range of 1.5 x 10~^ < Mgat/Mhaio < 0.02 is consistent with 
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approximately 100 satellites according to measured numbers in cosmological simulations. 
However, if we extend the mass range to lower and higher values there would be many more 
dark matter satellites, especially at the lower mass end. These satellites might potentially be 
detectable, add more stars to the tidal debris stellar halo, and thus change our conclusions. 
We therefore estimate the effects of the satellite mass range on the final total number of 
satellites that might be detected near M31. We sample the same mass function as equation 
(3) but this time with lower and upper limits 1.5 x 10~^ < Mgat/Mhaio < 0.1 that leads 
to about 1,000 satellites with a total combined mass of 0.25 Mhaio compared with a total 
satellite mass of 0.1 Mhaio for our 100 sateUite run. By extending the mass range, we add 
nearly 900 low mass halos and five high mass ones outside the original range. The six most 
massive halos would account for 0.16 Mhaio and would be easily detected having a central 
surface brightness greater than 27 mag arcsec"^ even with significant stripping as seen in 
our simulation. However, the 900 additional low mass halos would probably be fainter than 
32 mag arcsec^^, the central surface brightness of the lowest mass satellites in our run. More 
than 60% of them have a mass less than 2.5 x 10^^ and have central brightnesses of 34 mag 
arcsec"^ according to a linear extrapolation of our results. None of this population would 
be detected in the current surveys of the M31 stellar halo region so our exclusion of the low 
mass tail of satellite population has little effect on the predictions for the number of satellites 
aroTind M31. However, the low mass satellites contain a significant amount of mass and tidal 
stripping of this population could approximately double the amount of tidal debris in the 
stellar halo. Surface brightness maps might then be expected to be a magnitude higher than 
what our 45M particle simulation predicts. 



5. Summary &: Conclusion 

This study presents the evolution of a self-consistent population of satellites in the 
presence of a parent galaxy. The internal structure, mass function and spatial and orbital 
distribution of the satellite distribution is motivated by cosmological simulations. The work 
improves upon previous self-consistent studies (e.g. Font et al. 2001) by using a more realistic 
galaxy model and much improved numerical resolution. It follows the lead of current work on 
satellite tidal disruption to make quantitative predictions of the tidal debris field of galaxies 
(e.g. Johnston et al. 2001). 

This work shows that members of a typical population of subhalos orbiting in a galactic 
will come close to the disk with a ten billion years timescale. In addition, it is possible to 
model these satellites and make direct predictions on the number one would expect to detect 
around a typical galaxy. We also agree that it is possible to maintain a disk in spite of 
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substructure having a mass fraction of about O.lMhaio- 
Here we summarize the main results of this paper. 

1. Dynamical friction plays only a minor role in the evolution of the satellite system and 
the number density profile of the system is relatively unchanged over 10 Gyr. The 
orbits of the most massive satellites do show some decay but since they suffer severe 
tidal stripping, dynamical friction quickly becomes unimportant. 

2. The vast majority of the satellites survive in the galactic environment for more than 
a Hubble time; only 7 are completely disrupted. However, most of sateUites lose a 
significant fraction of their mass due to tidal interactions. 

3. The satellite mass loss due to tidal stripping is described by two distinct phases. The 
first one is characterized by a sharp mass decline with ti/2 = 3.5 Gyr. During that 
period, a small amount of mass loss, 5.3 x 10^ Mq, is due to a transient state where 
satellites initially located inside 50 kpc are overfilling their tidal radius and start losing 
mass instantly. During the second phase, the satellites lose about 25% of their initial 
mass with ti/2 = 9.3 Gyr. Over 10 Gyr, the satellite population mass diminishes by 
about 75% and turns into tidal debris. These debris form long tidal tails around the 
galaxy and can be detected in photometric maps of M31. The final density profile of 
the light satellites can be approximated by the fitting function given in Hayashi et al. 
(2003). 

4. The spatial distribution of stars associated with infalling satellites follow a power law 
p oc r~^'^ at large radii. This result is in agreement with recent numerical studies and 
observations of the stellar halo in edge-on disk galaxies. 

5. No obvious Holmberg Effect is observed. 

6. The mock B-band photometric maps, computed under the assumption that the most 
tightly bound particles are kinematic tracers of the stars, show conspicuous features of 
bridges and tails comparable to what is actually observed. A star formation efficiency 
of 10% is necessary to match the morphological and photometric properties of our 
Giant Stream proxy with the real one. 

7. The number of satellites detected around M31 below fiy ~ 26.8 mag arcsec"^ can 
match our results with a star formation efficiency of 10%. For a value of 1%, the 
number of sateUites detected goes down to 0. A value of 5% can not be discarded since 
dissipation is not taken into account. Many of the satellites between 26.8 and 29.4 mag 
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arcsec ^ in a 5% star formation efficiency scenario would probably be detected under 
the limit of 26.8 mag arcsec"^ if dissipation was included in the simulations. 

8. A star formation efficiency of about 5% is necessary to account for the actual M31 
value of the surface brightness of the stellar halo measured at 55 kpc. Comparisons for 
radii larger than 55 kpc are observationally challenging because of the very low surface 
brightness in the outer parts of the galaxy. The simulated maps show a flattening in 

the inner part of the surface brightness profile (for R < 55 kpc, /is ~ 28 mag arcsec"^) 
showing that dynamical friction is unable to bring the satellites close to the center. 
Most of the mass is lost at larger radii. 

9. The formation of a bar around 5 Gyr is the result of interactions of satellites with the 
disk. This intriguing result will be discussed in a forthcoming paper. 
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Table 1. 


Andromeda (M31) model 


Component 


Parameter 


Value 


Disk 








Mass 


7.77 X IQio Mq 




Rd 


5.57 kpc 




ZD 


0.3 kpc 




Rout 


30 kpc 


Bulge 








Mass 


2, 88 X 10i° Mq 




n 


1.82 kpc 






460 km s-i 




cut-off 


0.929 


Halo 








Mass 


5.83 X 10" Mq 




rs 


12.93 kpc 




Vs 


337 km s-i 




c 


13.37 




edge radius 


232.3 kpc 




cut-off 


0.75 
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Table 2. Physical parameters of the sateUite models. 



Name 




Vs [kpc] 


Vs [km s ^] 


cutoflF 


M(< rtid/Msat) 


c 


NFWl 


1.5 X 10-4 


0.7 


50 


0.4 


0.97 


29.4 


NFW2 


2.5 X 10-4 


0.85 


56 


0.4 


0.96 


27.7 


NFW3 


3.5 X 10"* 


0.98 


62 


0.4 


0.92 


26.8 


NFW4 


4.5 X 10-* 


1.1 


65 


0.4 


0.91 


25.4 


NFW5 


5.5 X 10--* 


1.2 


70 


0.4 


0.95 


25.2 


NFW6 


6.5 X 10--* 


1.33 


75 


0.4 


0.94 


24.5 


NFW7 


7.5 X lO""" 


1.39 


77 


0.4 


0.94 


24.2 


NFW8 


8.5 X 10"* 


1.46 


80 


0.4 


0.95 


24.0 


NFW9 


9.5 X 10-4 


1.54 


82 


0.4 


0.94 


23.5 


NFWIO 


1.5 X 10-3 


1.79 


93 


0.4 


0.93 


23.0 


NFWll 


2.5 X 10-3 


2.2 


112 


0.4 


0.91 


22.7 


NFW12 


3.5 X 10-3 


2.5 


122 


0.4 


0.89 


21.9 


NFW13 


5.5 X 10-3 


3.09 


141 


0.4 


0.88 


20.8 


NFW14 


6.5 X 10-3 


3.35 


146 


0.4 


0.87 


20.1 


NFW15 


2 X 10-2 


5.5 


200 


0.4 


0.77 


17.5 
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Table 3. B-band photometric maps parameters 



Kinematic tracers : 

Baryonic mass fraction (Mi^/Mdm) '■ 
Star formation efficiency'^ (Af*/Afb) '■ 
M31 M/Lb (constant)'' : 7.6 

For a simple stellar population formed at t = -1'' Gyr : 
M/Lb ratio at t=3.5 Gyr : 
M/Lb ratio at t=5.5 Gyr : 
M/Lb ratio at t=9.5 Gyr : 

For a simple stellar population formed at t = -2.5 Gyr : 
M/Lb ratio at t=3.5 Gyr 
M/Lb ratio at t=5.5 Gyr : 
M/Lb ratio at t=9.5 Gyr 
For a simple stellar population formed at t = -3.5 Gyr 
M/Lb ratio at t=3.5 Gyr : 
M/Lb ratio at t=5.5 Gyr : 
M/Lb ratio at t=9.5 Gyr : 



10 % most bound particles 
0.171" 
0.01, 0.05 and 0.1 



1.5476<^ 
2.0771 
2.9688 

1.9485 
2.4115 
3.2855 

2.1877 
2.6631 
3.4411 



''Steidel et al. (2003) 

''Ricotti & Gnedin (2005) 

■^Faber & Gallagher (1979) 

^GALAXEV code, Bruzual & Chariot (2003) 

°Notc that the simulation starts at t=0. We imply here that the stars formed prior to 
the beginning of the simulation. 
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Table 4. Number of satellite galaxies detected in a 10° x 10° field centered on M31 for 
different stellar population ages, star formation efficiency and isophotal thresholds 



tform (Gyr) 


Ceas (Gyr) 


M,/Mbar 


TV ^ 26.89 mag arcsec"^ 


N ^ 29.4 mag arcsec 


N ^ 31.9 mag 




3.5 


0.01 


1 


23 


26 




3.5 


0.05 


11 


26 


26 




3.5 


0.1 


25 


26 


26 




5.5 


0.01 





15 


30 




5.5 


0.05 


6 


30 


30 




5.5 


0.1 


15 


30 


30 




9.5 


0.01 





8 


29 




9.5 


0.05 


4 


30 


30 




9.5 


0.1 


8 


30 


30 


-2.5 


3.5 


0.01 





23 


26 


-2.5 


3.5 


0.05 


10 


26 


26 


-2.5 


3.5 


0.1 


23 


26 


26 


-2.5 


5.5 


0.01 





13 


30 


-2.5 


5.5 


0.05 


6 


30 


30 


-2.5 


5.5 


0.1 


15 


30 


30 


-2.5 


9.5 


0.01 





8 


29 


-2.5 


9.5 


0.05 


3 


30 


30 


-2.5 


9.5 


0.1 


7 


30 


30 


-3.5 


3.5 


0.01 


1 


20 


26 


-3.5 


3.5 


0.05 


9 
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The time at which the measurements are made after the beginning of the simulation. 
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Fig. 1. — Evolution of the disk velocity dispersion for the 35M control simulation. From 
bottom to top : a^. (J4, + 50 km s~^ and + 100 km s~^. Note that there is no vertical 
heating over 10 Gyr. Most of the heating in the radial and azimuthal directions occurs in 
the first 5 Gyr and is due to transient spiral features present early in the simulation. 
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Fig. 2. — Disk scale height versus radius for the control experiment. Scale height is taken to 
be the variance of z averaged over circular rings in the disk. Note that there is only minor 
growth of the scale height with time reflecting the quiet equilibrium of the disk DF and 
adequate resolution of the simulation. 
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Fig. 3. — Satellite cumulative mass function (top) and mass function (bottom). For the top 
figure, the plus symbols are the data points for our mass function and the filled line is the 
best dn/dMsat oc M~J^'^ fit to the data. Note that the mass function we use is normalized in 
such a way that in the interval 1.5 x 10~^ < Mgat/Mhaio < 0.02 the total mass in subhalos 
is 0.1 Mhaio- The other lines with symbols are taken from Gao et al. (2004). The X symbols 
and short-dashed line represent the subhalo mass function of a high resolution 2 x 10^^ 
Mq halo (GA3n run). The star symbols and the dotted line refer to an average of all the 
cluster- type simulations. 
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Fig. 4. — Evolution of the satellite cumulative number density function using solid test 
particles in the rigid M31 potential (no dynamical friction/tidal stripping). The different 
lines with symbols are for : plus-filled Gyr, X-dashed 5 Gyr, and stars-dotted 10 Gyr 
snapshot. The cumulative number density function evolves slightly over 10 Gyr reflecting the 
small number statistics and the fact that spherical NFW potential is only an approximation 
to the full potential for the galactic model. 



-32- 



20 



15 - 



10 - 



5 - 











50 100 150 200 250 300 

r [kpc] 



Fig. 5. — Distribution of the initial pericenter and apocenter passages radii for the sateUites 
population. Filled line : the pericenter radii; dashed line : apocenter radii. The pericenter 
distribution peaks at 50 kpc and the apocenter one at 250 kpc. The ratio ra/vp ~ 5 is typical 
of cosmological simulations (e.g. Moore et al. (1999)). 



-33- 





ii 










1 I 1 

■ 






L 






s 




i 




■ 


■ V 


- -P 

■ 


r 
m 


1 


» ■ 


i \ 


if 




ilk 







Fig. 6. — Snapshots of the 100-sateUite simulation. Prom top to bottom, t=0, 2, 4, 6, 8, and 
10 Gyr. One can see that a strong bar is formed between 4 and 6 Gyr. Several conspicuous 
shell structures arc visible, especially in the first few billion years. ((See high-resolution color 
snapshots at http://www.cita.utoronto.ca/~jgauthier/m31) 
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Fig. 7. — Evolution of the disk velocity dispersion for the simulation with subhalos. From 
bottom to top : az, + km and cr^ + 100 km s~^. Between 4 and 5 Gyr, a bar forms. 
This explains the sudden increase of ar and in the inner region of the galaxy between 3 
and 5 Gyr. 
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Fig. 8. — Same as Figure 2 but with satellite population. The scale height only grows slightly 
during the first 4 biUion years reflecting only mild heating from the sateUite population. The 
rapid heating thereafter is due to the formation of the bar. 
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Fig. 9. — Decay and mass loss of individual satellites. The first two columns are for the five 
most massive satellites. The third and fourth ones are for five of the lightest satellites. First 
and third columns give the radius of the center of mass as a function of time. The figure 
demonstrates that dynamical friction does not play much of a role in evolution of the orbits. 
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Fig. 10. — Evolution of the number density profile of satellites. The slope of the cusp remains 
constant over 10 Gyr.. 
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Fig. 11. — Density profiles for three different satellite mass bins at 10 Gyr. Top : 

Msat/Mhaio = 1-5 X 10-^ middle : Msat/Mhaio = 5.5 x 10"^ and bottom M^at/Mhaio = 
1.5 X 10~^. For the left plots, the thick line represents the unstripped density profile at the 
beginning of the simulation and the dashed line is the pure NFW model used to model the 
satellite. Arrows indicate the position of the scale radius r^. Satellites are initially truncated 
so that Tedge ^ ?^tidai- The right plots show the fit given by equation (8) for a typical profile 
of one of our sateUites in each mass bin. The short dashed line is the best fit and the dotted 
fine is the profile. As for the plots on the left, the dashed fine is the pure NFW model. 



-39- 




01 23456789 10 

t [Gyr] 

Fig. 12. — Evolution of the total mass bound in satellites as a function of time. The short- 
dashed and dotted lines are for the 45M and 450K simulation respectively. Long-dashed and 
filled lines are exponential decay fits to the 45M-particle run. The overall satellite decay is 
characterized by two distinct phases. The first one spans over the first 4 Gyr and corresponds 
to a very sharp mass loss with ti/2 = 3.5 Gyr. For the second one, from 4 to 10 Gyr, ti/2 ~ 
9.3 Gyr. 
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Fig. 13. — Number of satellites as a function of cos(^). Large fluctuations are initially present 
because of Poisson noise (y^lO) = 3) 
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Fig. 14. — 3D density profile of the stripped stars plus stars still in satellites after 10 
Gyr. We assume a constant star formation efficiency M^/Mb=0.1 and a baryon fraction of 
0.171. The long-dashed line represents the dark halo density profile; short-dashed one the 
stellar component of the bulge+disk; the dotted fine is associated with sateUites star only 
and the dash-dotted line is the sum of the disk, bulge and satellites stellar component. For 
r/r2oo > 0.2, the profile of sateUites stars is well fitted by p oc r~^-^. The spikes in the satellite 
stellar component are associated with stars that are still bound to orbiting satellites. 
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Fig. 15. — Scatter plot of the final mass nif vs initial mass rrii fo the 100 satellites in the 
experiment. Different symbols distinguish satellites according to pericenter. The arrows 
represent the seven satellites that are completely destroyed by the end of the simulation. 
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Fig. 16. — B-band surface brightness profiles of the stars tidally stripped from the satellites. 
The profiles are drawn assuming a star formation efficiency of 10% and a simple stellar 
population formed 1 Gyr prior to the beginning of the simulation. On both plots, the 
symbols are for four different snapshots : plus 2.5 Gyr, X 5.5 Gyr, stars 7.5 Gyr and squares 
10 Gyr. On the left-hand side, the dot-dashed line shows the contribution of disk+bulge of 
M31 to the total surface brightness. On the right-hand side, a closeup shows straight hnes 
corresponding to the best fit de Vaucouleurs profiles at 5, 7.5, and 10 Gyr. The feature 
around log R = 2.2 for the 2.5-Gyr snapshot is a consequence of the transient start and 
tends to disappear before 5 Gyr. 
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Fig. 17. — B-band CCD Photometric maps at 2000 x 2000 pixel resolution showing a 10° x 
10° field of the satellite stellar population taken at from top to bottom 3.5, 5.5 and 9.5 Gyr. 
For each map, we assume a simple stellar population formed at t=-2.5 Gyr. We show the 
contours corresponding to a surface brightness of 34.4, 31.9, 26.89, and 24.59 mag arcsec"^. 
From one figure to another, we vary the star formation efficiency i.e. how much gas is 
converted into stars, left column : Mgtars/^baryons = 0.01, middle : Mgtars/^baryons = 0.05 
and right column : Mstars/Msaryons = 0.1. We assume a constant Mbaryons/-^sat = 0.171 for 
all sateUites. The white patch in the center of M31 is to avoid saturation of the chip and 
corresponds to a limiting isophote of 24.59 mag arcsec"^. {See high-resolution color version 
of these maps at http://www.cita.utoronto.ca/~jgauthier/m31) 
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Fig. 18. — Close-up version of two Figure 17 maps showing features resembling both mor- 
phologically and photometrically to the Giant Stream and Gl Clump observed around M31. 
The surface brightness of our Giant Stream is about 28.5 mag arcsec"^ in the B-band.(S'ee 
high-resolution color version of these maps at http://www.cita.utoronto.ca/~jgauthier/m31) 



